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ABSTRACT 

We report on the observation of anisotropy in the arrival direction distribution of cosmic rays at 
PeV energies. The analysis is based on data taken between 2009 and 2012 with the IceTop air shower 
array at the South Pole. IceTop, an integral part of the IceCube detector, is sensitive to cosmic rays 
between 100 TeV and 1 EeV. With the current size of the IceTop data set, searches for anisotropy at 
the 10 -3 level can, for the first time, be extended to PeV energies. We divide the data set into two 
parts with median energies of 400 TeV and 2 PeV, respectively. In the low energy band, we observe a 
strong deficit with an angular size of about 30° and an amplitude of (— 1.58±0.46 s t a t ±0.52 sys ) x 10 -3 
at a location consistent with previous observations of cosmic rays with the IceCube neutrino detector. 
The study of the high energy band shows that the anisotropy persists to PeV energies and increases 
in amplitude to (-3.11 ± 0.38 sta t ± 0.96 sys ) x 1(T 3 . 
Subject headings: astroparticle physics — cosmic rays 
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1. INTRODUCTION AND MOTIVATION 

Over the last half-decade, several experiments in 
the northern and southern hemisphere have reported 
anisotropy in the arrival direction distribution of cosmic 
rays at TeV energies. In the northern sky, two features 
dominate the TeV cosmic ray sky: a large-scale struc- 
ture with a n amplitude of about 10 ~ 3 usually described 



as a dipole (Munakata et al. 1997 Amenomori et al. 2005 
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2006; |Guillian et al |2007||Abdo et aLp 009), and a small- 
scale structufe^itfT size 10° 
to 30° flAbdo et al |2008| |Vernetto et al.|2009[ ). 

In the southern hemisphere, cosmic ray arrival direc- 
tions at TeV energies observed with the IceCube neu- 
trino detector at the South Pole exhibit features simi- 
lar to those discovered in the northern sky. In a data 
set with a median energy of 20 TeV, the arrival direction 
distribution exhibits a large-scale structure similar in ori- 
entation and shap e to the large-scale fe ature observed in 
the northern sky ( |Abbasi et al.|2010b| ). 

In addition, there is a small-scale structure which is 
about a factor of five weaker in relative intensity than the 
large-scale structure. This small-scale structure contains 
several regions of s ignificant cosmic-ray excess and deficit 
( |Abbasi et aLl|2011[ ). 

There are several models that can at least qualitatively 
explain the anisotropy. Cosmic rays in this energy range 
are assumed to be accelerated in Galactic sources, most 
likely in shocks from supernova explosions. The trans- 
port of cosmic rays at these energies in the Galactic mag- 
netic field is diffusive, and the flux from a single nearby 
source would be observed on Earth as a dipole with its 
maximum possibly oriented towards the source. If a few 
supernova remnants from recent (10... 100 kyrs) nearby 
supernovae were primarily responsible for t he Galactic 
cosmic ray flux (Erl ykin fc Wolfendale|2006l ), their com- 
bined flux on Earth would be a superposition of these 
individual dipoles. The observed large-scale structure 
in the cosmic ray flux could be the sum of the contri- 
butions from a few nearby sources and from the large 
scale distribution of supernova remnants in our Galaxy 



(Blasi & Amato 2012 Pohl & Eichler 2012). Our lim 



ited knowledge oil nearby supernova remnants renders a 
more quantitative explanation of the amplitude and the 
phase of the observed large-scale anisotropy impossible, 
but among the predictions of this model is an increase 
of the anisotropy with the energy of the primary cosmic 
rays. 

The small-scale structure is more difficult to explain. 
The Larmor radius of a proton with 10 TeV energy in 
a magnetic field with jiG strength (an estim ate of the 
strength of magnetic fields in our Galaxy (Han et al. 
2006)) is only of order 0.01 pc, so the observed small- 
scale anisotropy cannot correlate with nearby sources. 
Recent studies link the smaller structure to cosmic ray 
progation in turbulent magnetic fields within a few tens 
of parsecs from Earth (Giacinti & Sigl 2012). In this 
case, the small-scale structure is also expected to show a 
dependence on energy. 

Several experiments have studied the energy depen- 
dence of the anisotropy. In the northern hemisphere, the 
EAS-TOP experiment reports anisotropy up to at least 
400 TeV. The data contain weak evidence for an increase 
in the amplitude of the anisotr opy as a function of energy, 
as well as a change of phase (Agli etta et al.||2009 ). Mea- 
surements in the southern hemisphere with IceCube also 
indicate that anisotropy is still present in a data set with 
400 TeV median energy. It differs in shape and strength 
from the anisotropy observed at 20 TeV and is no longer 
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a superposition of large- and small-scale structures, but 
rather dominated by a single deficit region with an an- 
gular size of about 30° flAbbasi et aL|2012b| ). 

IceCube is primarily a neutrino detector designed to 
search for sources of astrophysical neutrinos. The data 
set used in the cosmic ray analysis consists of downgoing 
atmospheric muons from cosmic ray air showers in the 
atmosphere above the IceCube detector. The downgo- 
ing muons preserve the direction of the cosmic ray pri- 
mary, but the muon energy is only a poor indicator of 
the air shower energy. In addition, because of the high 
muon rate (10 6 times the neutrino rate in IceCube), these 
events are stored in a separate data format which only 
contains the results of a fast online reconstruction per- 
formed at the South Pole. No raw data are preserved. 
An offline reconstruction with more sophisticated algo- 
rithms to determine the muon energy is therefore not 
possible. 

Abov e 1 TeV, the energy r esolution is of order 0.5 in 
\og(E) ( |Abbasi et aL||2012b| >- The energy resolution is 
estimated as the standard deviation of the distribution 
of the difference log(E true ) — log(£? reco ), where E true and 
E reco are the true and reconstructed shower energies ob- 
tained from simulation studies, respectively. This distri- 
bution has substantial tails, making it difficult to isolate 
a set of events with large median energy that is not con- 
taminated by low-energy events. 

The IceTop air shower array is located at 2835 m alti- 
tude on the surface of the ice sheet above the IceCube 
neutrino detector. IceTop is a dedicated cosmic ray de- 
tector optimized for air shower observations at PeV ener- 
gies. IceTop is used to record not only the muonic com- 
ponent of the air showers, but also the electromagnetic 
component, at ground level. With the sparse sampling 
of the shower front typical for air shower arrays, it also 
has a considerably higher detection threshold for cosmic 
rays than IceCube. The size and geometry of the array 
result in a threshold for reconstruction of air showers of 
approximately 300 TeV. 

As an air shower array, IceTop provides a more mea- 
sured information per shower than IceCube. A study of 
cosmic ray anisotropy with IceTop can therefore com- 
plement measurements with IceCube. With its high en- 
ergy th reshold, an energy re solution better than 0.1 in 



\og(E) (Abbasi et al. [2012a ), and sensitivity to the cos- 
mic ray composition, IceTop data are particularly useful 
for studying anisotropy at energies above 1 PeV. Due to 
the lower data rate, the full event information is stored, 
and a more careful offline reconstruction of the primary 
cosmic ray properties is possible. 

The accumulated IceTop data set is not yet large 
enough for a detailed study of per-mille anisotropy in 
several energy bins. However, the statistics collected 
between 2009 and 2012 are now sufficient to search for 
anisotropy in two energy bands centered at 400 TeV and 
2 PeV. The 400 TeV data set can be compared to results 
based on downgoing muons in IceCube at a similar en- 
ergy but for a differe nt cosmic ray composition model 
flAbbasi et aT]|2012b| ). With the 2 PeV data set, the 
search lor anisotropy is extended to energies not pre- 
viously explored. In this paper, we report on the obser- 
vation of anisotropy with IceTop at both energies. 

The paper is organized as follows. In Section 2, we de- 
scribe the IceTop detector and the data sample used in 



this analysis. The analysis techniques and the results are 
briefly described in Section 3. The analysis is based on 
methods used in previous work; a more de tailed descrip- 
tion o f these techniques can be found in | Abbasi et al. 
fl2011[ ) . Section 4 discusses systematic uncertainties, and 
Section 5 summarizes the paper. 

2. DETECTOR, DATA SETS AND SIMULATION 
2.1. The IceTop Detector 

The IceTop cosmic ray air shower array consists of 81 
stations distributed over an area of 1 km 2 in a hexagonal 
grid with a distance of about 125 m between neighbor- 
ing stations. During the construction phase of IceTop 
between 2005 and 2010, the detector was operated in 
several partial configurations. In this work we use data 
taken during three periods: between May 2009 and May 
2010 when the detector was operated in a 59-station con- 
figuration (IT59); between May 2010 and May 2011 when 
IceTop operated with 73 stations (IT73); and between 
May 2011 and May 2012 when the detector operated in 
its final 81-station configuration (IT81). The layout of 
these configurations is shown in Fig. [T] 

IceTop is described in detail in |Ab"r5asi et al.| fl2012a[ ) . 
Each IceTop station is instrumented with two light-tight 
ice Cherenkov tanks separated by about 10 m. Each tank 
is 1.8 m in diameter, 1.3 m in height, and is filled with 
transparent ice up to a height of 0.9 m. Frozen into 
the ice are two Digital Optical Modules (DOMs) that 
are used to detect the Cherenkov radiation emitted by 
charged leptons present in the cosmic ray air shower. 
A DOM consists of a glass sphere that houses a 10" 
Hamamatsu photomultiplier tube (PMT) ( [Abbasi et al. 
|2010a| ) , together with elec tronic boards used l or filtering, 
digitization, and readout ( Abbasi et al.|[2009 ). 

The two DOMs inside each IceTop tank are operated 
at different PMT gains in order to increase the dynamic 
range of the detector. The high-gain DOMs in the two 
tanks that form a station are run in local coincidence 
mode, and data readout is enabled if one of the DOMs 
records photon hits within ±1 /is of the other. The trig- 
ger condition in IceTop requires at least six DOMs to 
have recorded locally-coincident hits within a time win- 
dow of 5 /is, which implies that at least two stations 
participated in the event. 

Due to the limited bandwidth available for data trans- 
mission from the South Pole, events triggering less than 
eight stations were prescaled by a factor of eight during 
the operation of IT59, and by a factor of three during 
IT73 and IT81. Events triggering more than eight sta- 
tions were not prescaled. 

2.2. Data Sample and Simulation 

The prescaling scheme described above was used to 
divide the data into two samples: a "low-energy" data 
set, containing events with at least three but less than 
eight stations triggered, and a "high-energy" data set 
that contains events where eight or more stations were 
triggered. 

During the operation of IT59, IT73, and the first year 
of IT81, a total of 3.55 x 10 8 events with more than 3 trig- 
gered stations were recorded. Of these events, 2.90 x 10 8 
were classified as "low-energy" events while the "high- 
energy" sample contains the remaining 0.65 x 10 8 events. 
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A zenith angle cut (described below) was used to remove 
misreconstructed events at large zenith angles. This cut 
reduced the final sample to 2.86 x 10 8 events in the low- 
energy set and 0.64 x 10 8 events in the high-energy sam- 
ple. 

The angular resolution of the shower reconstruction al- 
gorithm and the median energy of the data sets were de- 
termined using simulated cosmic ray air showers. Events 
were generated w ith the CORSIKA Monte Carlo code 
(Hec k et "aLp 998) and passed through a f ull simulation 
ol the IceTop detector ( |Abbasi et al.|2012a| . The median 
energy of the samples determined using this simulation 
will depend on the assumptions made about the chemi- 
cal composition of the primary cosmic rays. The detailed 
primary composition has not been directly measured for 
energies beyond 100 TeV, but models that extrapolate 
existing measurements to higher energies indicate that 
in the energy range of this analysis, the cosmic ray flux 
consists mainly of protons, helium, and iron (Gaisser 



2011 ). Their relative contribution is a function of energy, 
with helium and protons dominating around 100 TeV and 
iron becoming the dominant element above several tens 
of PeV. Given the uncertainties in the composition, we 
have generated only proton and iron showers as the two 
limiting cases for the chemical composition. The true 
median energy of the sample should be contained in the 
interval defined by these two cases. 

The arrival direction of cosmic ray showers in IceTop 
is reconstructed by fitting a plane to the front of the air 
shower. The fit algorithm implements an analytic x 2 - 
minimization that uses the positions and hit times of the 
triggered stations to reconstruct the direction vector of 
the shower. From simulation we have determined the 
median angular resolution of this algorithm to be 3° for 
both proton and iron showers for all detector configura- 
tions. In other IceTop analyses, a more precise recon- 
struction algorithm is used that takes into account the 
curvature of the air shower front and can reach a sub- 
degree resolution. The plane fit is better suited to our 
needs since it provides a resolution that is several times 
smaller than the typical angular scale of the anisotropic 
pattern (> 20°) without requiring a larger number of sta- 
tions triggered which would reduce the size of the cosmic 
ray sample. As shown in Fig. |2j the resolution of the 
plane fit degrades rapidly for showers with zenith angles 
larger than 60°. For this reason, this analysis is limited 
to events with a reconstructed zenith angle smaller than 
55°. 

The median energies of the data sets were determined 
from the energy distribution of the simulated air show- 
ers which satisfy the same trigger conditions as events in 
the low- and high-energy data samples. The simulated 
energy distributions are shown in Fig(3j The median 
energies and the 68% containing intervals for the two 
composition models for each energy band are shown in 
Table [I] which shows that the low-energy band has a me- 
dian energy in the range 270 - 500 TeV, while the median 
energy for the high-energy sample should be contained in 
the range 1.6 - 2.2 PeV. 

For both data samples, the median energy of the pri- 
mary cosmic rays monotonically increases with zenith- 
angle, as illustrated in Fig. [4j 

3. DATA ANALYSIS AND RESULTS 
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Figure 1. Detector configurations of the IceTop array, 2009-2011. 
IT59 comprised 59 stations deployed through January 2009 (blue 
circles). In 2009 and 2010 fourteen additional stations were de- 
ployed and the detector was operated in the IT73 configuration 
(blue and green circles). The remaining eight stations (red circles) 
were deployed in late 2010. The final IceTop configuration (IT81) 
consists of 81 stations and operated in 2011. 
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Figure 2. Median opening angle between reconstructed and 
true arrival direction as a function of reconstructed zenith angle 
0. At large zenith angles the fraction of misreconstructed events 
increases. For this reason, a zenith cut was implemented that re- 
stricts the analysis to events with < 55° (dashed gray line). The 
error bars correspond to a 68% containing interval. IT59, IT73, 
and IT81 show the same dependence of angular resolution on re- 
constructed zenith angle. 



3.1. 



Map Making Procedure 

The search for anisotropy in the IceTop data is based 
on techniques applied in the analysis of cosmic ray 
with IceCube and described in more detail in 



r ay data 
lAbb asi 



et al. (2011 ). To obtain a skymap of the relative intensity 
ol the cosmic ray flux, the map of reconstructed arrival 
directions in equatorial coordinates is compared to an 
estimate of the isotropic expectation represented by a 
"reference map." The reference map is cons tructed from 
the data us ing the time-scrambling method ( |Alexandreas| 
et al. 1993[ ). For each detected event, 20 "fake" events are 
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Low 


energy 


High energy 


Composition 


E 


68% interval 


E 68% interval 


Proton 


0.27 PeV 


0.11- 0.69 PeV 


1.6 PeV 0.83 - 3.8 PeV 


Iron 


0.50 PeV 


0.22 - 1.2 PeV 


2.2 PeV 1.2 - 5.3 PeV 



Table 1 

Median energy and 68% containing interval in PeV for the two energy bands used in this work assuming that the cosmic rays consists of 

either protons or iron nuclei. 
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Figure 3. Simulated energy distributions for all events in the low- 
energy {dashed) and high-energy (solid) data sets assuming all-iron 
(blue) and all-proton (red) compositions. The energy distributions 
are the same for the IT59, IT73 and IT81 configurations. 

generated by keeping the local zenith and azimuth angles 
fixed and calculating new values for equatorial coordi- 
nates using event times randomly selected from within a 
time window At bracketing the event. The fake events 
are stored in the reference map with a weight of 1/20. In 
order to be sensitive to anisotropy at all angular scales, 
we use At = 24 hr. The stability of the detector over 
this time was verified by performing a x 2 -test where the 
distributions of zenith and azimuth coordinates for the 
events are compared inside the window. 

By using the time-scrambling algorithm, the events in- 
cluded in the reference map have the same arrival di- 
rection distribution as the data in local coordinates. In 
addition, the method compensates for variations in the 
event rate, including gaps in the detector uptime. We 
note, however, that the method is known to create ar- 
tificial deficits near strong excesses, and vice versa near 
large deficits. With a 24 hour scrambling window, a sin- 
gle strong excess (deficit) will raise (lower) the reference 
map level for the entire right ascension band at the dec- 
lination of the excess (deficit) which could bias the ob- 
served amplitude of the anisotropy and its statistical sig- 
nificance. This effect can become important if extended 
regions of strong excess or deficit flux are present in the 
data, and the resulting skymaps should be interpreted 
carefully with this potential bias in mind. 

The maps are built u sing the HEALPix e qual- area pix- 
elisation of the sphere ( Gorski et al.||2QQ5 ) with an aver- 
age pixel size of about 0.9 U . Maps of statistical signifi- 
cance for the low- and high-energy data sets are shown 
m Fig. [5] The bin size is not optimized for a study of 
anisotropy at scales larger than the angular resolution of 
the detector. We therefore apply a smoothing procedure 
to both the data map and the reference map in order 
to increase the sensitivity of the method to structures 



with larger angular sizes. The smoothing procedure is 
essentially a rebinning of the maps, but rather than pro- 
ducing maps with fewer (but statistically independent) 
pixels, we retain the original 0.9° binning. At each bin, 
we add the counts from all pixels within some angular 
radius ("smoothing radius") of the bin. This produces 
maps with Poisson uncertainties, but with bins that are 
not statistically independent. Since the optimal smooth- 
ing scale is not known a priori, we perform a scan over 
different smoothing angles. After smoothing, the relative 
intensity 5Ii, i.e. the amplitude of deviations from the 
isotropic expectation for each angular bin z, is calculated 
using: 

AN, Ni - (N)i 



where Ni and (N)i are the number of events in pixel 
i of the data map and the reference map, respectively. 
The statistical significance of the obser ved devia t ions i s 
calculated using the method described in |Li fc Ma ( 1983). 

The significance is later corrected to account tor the 
number of trials introduced by looking everywhere in 
the sky for significant fluctuations, and for the fact that 
a scan was performed over different smoothing radii to 
search for anisotropy at different angular scales. 

3.2. Results 

The smoothing procedure described above was applied 
to the low- and high-energy maps for smoothing radii be- 
tween 5° and 50° in 3° steps. A search for regions of high 
significance was performed on the resulting smoothed 
maps. The relative intensity and significance maps for 
the low- and high-energy data are shown in Fig. 6] for 
a representative smoothing radius of 20° where all the 
relevant features observed in these two energy ranges are 
visible. Maps for all other smoothing radii are available 
on the web as supplemental information to this paper. 

The low-energy map is dominated by a strong deficit 
in relative intensity. The statistical significance of the 
deficit reaches a maximum of 8.5<r for a smoothing ra- 
dius of 29° at a location around (a = 85.8°, S = —36.4°). 
Since the search for this minimum is performed over 
about 10000 pixels in the map, and across all 16 dif- 
ferent smoothing radii, there is a trials factor of at most 
1.6 x 10 5 that reduces the post-trial significance of the 
deficit to 7.0<j. It must be noted that this correction for 
trials is conservative, since the pixels in the map are sta- 
tistically correlated by the smoothing procedure, which 
results in a smaller effective number of trials than the 
maximum. 

For the optimal smoothing radius of 29°, the relative 
intensity SI reaches a value of about —1.5 x 10 -3 at 
the location of the greatest deficit around (a = 83.7°, 
5 = —35.7°), near the edge of our exposure window. Dif- 
ferences in declination between the location of the max- 
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Figure 4. Median energy as a function of reconstructed zenith angle for the low-energy (left), and high-energy (right) data sets for proton 
and iron cosmic ray primaries. The error bars correspond to a 68% containing interval. IT59, IT73 and IT81 show the same energy 
dependence on the reconstructed zenith angle. 
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Figure 5. Maps of statistical significance for the low-energy (top) 
and high-energy (bottom) data sets with no smoothing applied. 

imum relative intensity and maximum significance are 
due to the fact that the statistical significance accounts 
for both signal strength and the declination-dependence 
of our statistics. This usually implies that the position of 
maximum significance is offset towards lower declination 
values where the statistics increase. 

Also visible in the low-energy map is a region of excess 
flux located around (a = 182.9°, 5 = -55.9°). The 
maximum pre-trial significance for this region is 5.3(7 for 
a smoothing angle of 26°. The significance falls below 
the 3cr threshold after accounting for trials. 

As mentioned above, in the presence of a strong deficit 
the time-scrambling algorithm can introduce an underes- 
timation of the isotropic reference level, which can pro- 
duce spurious excess regions in the parts of the sky sur- 
rounding the signal region. For this reason, it is possible 
that the excess observed in the low-energy data set is 
associated with the presence of the deficit region. 

The high-energy map also shows statistically signifi- 
cant anisotropy which is dominated by a deficit located 
in the same approximate position as that observed in the 
low-energy data. The pre-trial significance of the deficit 



is 8.6a (7.1a post-trials) for a smoothing angle of 35°, 
with its minimum located at (a = 79.4°, 5 = —37.2°). 
The main difference between the low- and high-energy 
deficits is that the value of 51 for the greatest deficit in 
the high-energy sample is —2.3 x 10 -3 , larger than its 
low-energy counterpart. This is evident in Fig. [7| where 
the relative intensity is projected onto the right ascension 
axis using the declination band —75° < 5 < —35°. 

A second notable feature in the high-energy map is 
a wide excess region that reaches a peak significance of 
5.9a (3.4a post-trials) for a smoothing angle of 41°. The 
excess does not appear to be concentrated in any partic- 
ular part of the sky, but distributed across a wide band 
in right ascension. This is visible in the one-dimensional 
projection shown in Fig. [7| where the relative intensity 
reaches a plateau above a > 170° which is offset from 
zero by about 10 -3 . As in the low-energy data set, such 
an excess could be associated with the presence of the 
observed deficit in the same declination band that intro- 
duces a bias in the reference-level estimation. 

In order to characterize the observed anisotropic pat- 
tern, we attempted to fit the relative intensity projections 
of the data shown in Fig. [7] with the first terms (dipole 
and quadrupole) of a harmonic series: 

2 

51(a) = J2A j cos[j(a-h)] + B . (2) 
j=i 

The harmonic fit parameters for the low-energy dataset 
are A x = (5.53 ± 0.91) x HT 4 , fa = -111.7° ± 9.7°, 
A 2 = (4.03 ± 0.80) x 10~ 4 , fa = 1-0° ± 7.9°, and B = 
(—0.47 ± 0.66) x 10 -4 . For all parameters the indicated 
uncertainties are statistical. The x 2 /dof for the fit is 
18.4/10. 

In the case of the high-energy dataset the fit parame- 
ters are A x = (1.43 ± 0.19) x 10" 3 , fa = -84.7° ± 8.0°, 
A 2 = (7.34 ± 1.9) x 10- 4 , fa = -14.7° ± 7.7°, and 
B = (-0.14 ± 1.35) x 10~ 4 . In this case, the x 2 /dof 
for the fit is 6.3/10. 

In the case of the low-energy dataset, the reduced % 2 of 
the fit is considerably larger than unity, indicating that 
for this data, as in the case of the previo us observation 
of anis otropy at 400 TeV with IceCube ( |Abbasi et al. 
|2012b[ ) , this choice of harmonic base functions does not 
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fit the data particularly well. For this reason, a new fit 
is performed using the following Gaussian function: 



61(a) Ae 



-(■ 



\/2c 



(3) 



where a is right ascension, A is the amplitude, a is the 
width, and a s is the right ascension of the center of 
the deficit. The parameter b represents an overall off- 
set from isotropy that can be introduced by the presence 
of a strong signal in the data. 

The results of these fits are shown in Table [2j and 
indicate that while the center point of the deficit for both 
data sets is consistently located at a s ~ 90°, both the 
amplitude and the width are larger in the high-energy 
sample, with both values increasing by about a factor of 
two with respect to the low-energy case. The location 
of the deficit in the right ascension projection (Fig. [7]) 
is consistent with its location in the skymap (Fig. pF, 
within statistical and systematic uncertainties, for both 
the low-energy and high-energy samples. Similarly, the 
amplitudes in relative intensity of the deficit agree well, 
when the overall offset b is taken into account. 

The systematic uncertainty associated with each fit 
value was obtained from the systematic uncertainty of 
the relative intensity projection, shown as shaded boxes 
in Fig. [7J The systematic uncertainty of the projection 
was conservatively estimated as the maximum amplitude 
of the relative intensity distribution for each energy data 
set when analyzed in the anti-sidereal time frame (see 
Section [I]). 

A previous search for anisotropy as a function of 
cosmic-ray primary energy was p erformed using muon 



data from the IceCube detector dAbbasi et al. 
In that analysis, cuts were applie< 



2012bp . 

to the data to select 



two sets of cosmic-ray events: one with a median pri- 
mary energy of 20 TeV, and a second one with a median 
energy of 400 TeV, similar to the low-energy IceTop sam- 
ple. The 400 TeV IceCube skymap shows a deficit region 
similar to the one observed in the IceTop low-energy sam- 
ple. At 20 TeV, IceCube observed a large-scale structure 
that is consistent with previous observations at th ese en- 
ergies ( |Abbasi et aI1|2010b[ )( [Abbasi et al.||2011[ ). The 
angular size and the orientation in the sky of the 20 TeV 
anisotropy is different from that observed at 400 TeV by 
IceCube and IceTop. Note that while the IceTop median 
energy was obtained by assuming two limiting cases of 
primary chemical composition (all-proton or all-iron pri- 
maries only), the IceCube median energy was obtained 
assuming that the co smic ray flux fol lows the polygonato 
composition models ( jHorandeI||2003 ) , which in principle 
could lead to some differences in the actual median en- 
ergy and the energy distribution of events in both sam- 
ples. 

The smoothing procedure described here was also used 
in the IceCube analysis, and the significance of the deficit 
was maximized for a smoothing radius of 29° , the same as 
the optimal smoothing angle for the low-energy IceTop 
data. Fig. [8] shows a comparison between the IceCube 
and IceTop results at 400 TeV. The location and ampli- 
tude of the deficits observed in both data sets agree given 
the current values of the statistical and systematic un- 
certainties associated with both measurements. 



A number of tests have been performed in order to 
quantify the systematic uncertainties associated with the 
observation of anisotropy in the IceTop data. 

In the first study, the anisotropy search was performed 
on three independent data subsamples, each containing 
events recorded during the operation of the three differ- 
ent detector configurations IT59, IT73, and IT81 con- 
sidered in this work. In this manner we can determine 
the possible systematic effect introduced by the changing 
geometry of the detector on the observed anisotropy. 

The results of this comparison are shown in Fig. [9j 
where the relative intensity as a function of right ascen- 
sion for the declination band —75° < 6 < —35° is dis- 
played for all three detector configurations and for the 
low- and high-energy samples separately. The anisotropy 
observed by all three configurations is consistent within 
statistical uncertainties. 

Another test was performed to evaluate the impact of 
the s easonal varia tion of the cosmic ray rate at the South 
Pole ( Til av et al.|2009[ ) . In this study, four different time 
periods were selected from the data: June through Au- 
gust, September through November, December through 
February, and March through May for each year of op- 
eration of the detector. These four data sets contain 
events taken with comparable detector geometries, but 
recorded during different phases of the seasonal varia- 
tion cycle. The results of this study indicate that the 
anisotropy observed in each of the four time periods is 
consistent within statistical uncertainties. 

Other possible seasonal effects on the anisotropy are 
also investigated. First, an analysis was performed to 
look for anisotropy in the so-called "solar" time frame, 
defined as having 365.25 (i.e. complete revolutions in 
the coordinate frame) per year. The motion of the Earth 
around the Sun should create a dipolar anisotropy in 
the solar frame with an amplitude of 4.7 x 10 -4 . No 
anisotropy was observed using IceTop data. However, 
simulations of the solar dipole assuming the IceTop ac- 
ceptance in local coordinates indicate that the current 
size of the data set is insufficient for a statistically sig- 
nificant observation. 

The second analysis consists of a search for anisotropy 
analysis in the "anti-sidereal" time frame, defined as hav- 
ing 364.25 days . No signal should be observed in this 
frame unless there exists a seasonal variation in the solar 
time frame that could affect the anisotropy in s i dereal 
time (period of 366.25 days). See |Abbasi et al.j ( |2Q1 1| ) 
for details. 

We performed the anti-sidereal analysis on the com- 
bined three-year data set and obtained both skymaps 
and one-dimensional relative intensity projections for the 
low- and high-energy bands. The skymaps produced 
for the anti-sidereal frame do not exhibit any significant 
anisotropy that could indicate a possible systematic bias 
in the sidereal frame. The systematic uncertainty of the 
sidereal anisotropy due to seasonal variations, shown in 
Fig. [7j is obtained from the relative intensity projections 
in the anti-sidereal frame. This uncertainty is conser- 
vatively estimated as the maximum departure from the 
reference level of the anti-sidereal right ascension distri- 
bution. 
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Low energy High energy 

A (-1.58 ± 0.46 ± 0.52) x lO" 3 (-3.11 ± 0.38 ± 0.96) x lO" 3 

a s 90.6° ± 6.8° ± 9.3° 88.1° ± 6.8° ± 11.1° 

cr 21.3° ±5.8° ±7.6° 43.1° ±7.3° ± 13.1° 

b (2.61 ± 0.64 ± 5.20) x lO" 4 (9.37 ± 1.96 ± 9.60) x lO" 4 

X 2 /dof 13.2/11 10.7/11 



Table 2 

Fit parameters obtained for both energy datasets for the Gaussian function given in Eq. [3] In all cases, the first quoted uncertainty is 

statistical while the second one corresponds to the systematics. 
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Figure 6. Relative intensity (top) and statistical significance (bottom) maps for the low-energy (left) and high-energy (right) data sets 
for a smoothing angle of 20°. 



A study of cosmic ray arrival directions with IceTop at 
two different median energies, 400 TeV and 2 PeV, shows 
significant anisotropy in both sets. The skymap is dom- 
inated by a single deficit region with an angular size of 
about 30° . The skymap at 400 TeV is similar to a skymap 
of comparabl e median energy o btained from cosmic rays 
in IceCube ( |Abbasi et al.|2012b ). IceTop data show that 
this anisotropy persists to 2 PeV. 

The anisotropy in the southern sky at 400 TeV and 
2 PeV is different in shape and amplitude from what 
is observed at 20 TeV. In the northern hemisphere, the 
EAS-TOP experiment has also recently found indications 
for an increasing amplitude and a change of phase be- 
tween 100 TeV and 400 TeV in a harmonic analysis in 
right a scension that consid ers the first and second har- 
monic ( |Aglietta et al.||2009| ). The IceTop anisotropy is 
not well- described by a sum of a dipole and a quadrupole 
moment, so the results cannot be directly compared. 
However, both northern and southern hemisphere data 
seem to show qualitatively similar trends. 

Although these results do not provide conclusive evi- 
dence for any particular model, they lend support to sce- 
narios where the large-scale anisotropy is a superposition 
of the flux from a few nearby sources. The sparse spatial 
distribution and the different ages of nearby supernova 
remnants are expected to lead to a bumpy structure in 
the amplitude and sudden chang es in the phase of th e 
anisotropy as a function of energy flBlasifc Amato|2012p . 



Unfortunately, this energy dependence is dominated by 
details such as the geometry of the Galaxy, the location, 
age and injection spectrum of the sources, and the en- 
ergy dependence of the cosmic ray diffusion coefficient. 
While the predicted strength of the amplitude has the 
correct order of magnitude, further quantitative predic- 
tions are not possible at this point. In addition, in their 
simplest form, these models predict a dipolar anisotropy, 
whereas in most cases, the observed anisotropy cannot 
be described as a simple dipole, which also means that 
"amplitude" and "phase" are not well-defined. 

It was recently pointed out that an existing dipolar 
flux in addition to cosmic ray propagation in turbulent 
magnetic fields close to Eart h can explain the ap pear- 

3ijIpQ 12|), For 
\Tth£ 



le re 



evant 



ance of small-scale structure (Giacinti & Sig 
cosmic rays with energies from TeV to PeV, 
distance scale is a few tens of parsecs, so the observed 
anisotropy at these energies is indicative of the turbu- 
lent Galactic magnetic field within this distance from 
Earth. The model predicts that the anisotropy is energy- 
dependent, but again, due to our poor knowledge of inter- 
stellar magnetic fields, it cannot provide more quantita- 
tive predictions that can be tested with data. A detailed 
measurement of the anisotropy might lead to a better 
understanding of these fields. 

The observation of cosmic ray anisotropy with IceTop 
opens up new possibilities for future studies that go be- 
yond mapping the arrival direction distribution as a func- 
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Figure 7. Relative intensity as a function of right ascension for the low-energy (left) and high-energy (right) data samples in the declination 
band —75° < 5 < —35°. The error bars are statistical while the colored boxes indicate the systematic uncertainty obtained from analyzing 
the same data in the anti-sidereal time frame (see Section^] for details). The result of a fit using the Gaussian function given in Eq. [3] to 
both energy bands are also shown. 

Some of the result s in this paper hav e been derived us- 
ing the HEALPix ( jGorski et al. |2QQ5| ) software libraries. 
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Figure 8. Comparison between the relative intensity projections 
for the IceTop low-energy sample (blue filled circles) and the Ice- 
Cube 40 TeV sample (black open circles) reported by |Abbasi et aD 
(2012bl. The location and amplitude of both deficits are consistent 
given the statistical and systematic uncertainties. The declination 
range for the IceCube plot is —75° < 5 < —25°, slightly different 
from the IceTop one. 

tion of energy. IceTop is designed to measure the energy 
spectrum and the chemical composition of the cosmic 
ray flux above several hundred TeV, and these capabil- 
ities allow for additional studies of the anisotropy. For 
one of the excess regions observed in the 10 TeV sky map, 
the Milagro experiment has reported a diffe rent energy 
spect rum than the isotropic cosmic ray flux ( |Abdo et al. 
2008). With data from IceTop, studies of the energy 
spectrum and composition of the cosmic ray flux in dis- 
tinct regions of the southern sky can be performed. 

IceTop is now in stable running mode in its complete 
configuration of 81 stations. In two years, the size of 
the cosmic ray data set available for anisotropy studies 
will be more than twice what was used in the analysis 
presented in this paper. Eventually, it will be possible 
to extend the analysis of cosmic ray anisotropy to higher 
energies. 
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